module mod_myelg
implicit none



ABSTRACT INTERFACE
  FUNCTION function_template(x) RESULT(res)
      double precision, INTENT(in) :: x
      double precision :: res
  END FUNCTION function_template
  END INTERFACE

  double precision :: H0, OmegaM, z0, zfrom, zto, dz, wa, w0, xi, ns, OmegaB, gamma, gammaK, D0, CH0

  double precision :: deltaH, h, keq

  double precision, allocatable, dimension(:) :: z_bins

  integer :: no_zbins, intN, curi, curj, curl


contains

  
  subroutine myelg_initialize

    integer:: i, N

    zfrom = 0.01
    zto = 0.2
    
    H0 = 100
    h = 0.71
    OmegaM = 0.3
    OmegaB = 0.04
    z0 = 0.8/1.4
    no_zbins = 6
    dz = (zto - zfrom)/(no_zbins)
    wa = 0.00
    w0 = -1.0
    ns = 0.96

    CH0 = 2997.9

    intN = 2000
    
    allocate(z_bins(no_zbins))
    do i = 1, no_zbins
       z_bins(i) = zfrom + (zto-zfrom)/(no_zbins)*i
    end do

    z_bins(1)=0.01d0         
    z_bins(2)=0.50d0
    z_bins(3)=0.71d0
    z_bins(4)=0.91d0
    z_bins(5)=1.19d0 
    z_bins(6)=2.5d0
    
    

    call myelg_calculate_deltaH()

    gamma = 0.55 + 0.05*(1.0 + w0 + 0.5*wa)
    !write(*,*) "gamma:" , gamma
    gammaK = OmegaM*h*exp( - OmegaB * (1 + sqrt(2.0*h)/OmegaM))
    keq = 0.073*gammaK
    D0 = myelg_D(1.d0)


  end subroutine myelg_initialize

function myelg_calculate_shear(i, j, l)

  integer, intent(in) :: i,j ,l
  double precision :: myelg_calculate_shear
  PROCEDURE(function_template), POINTER :: func

  func => myelg_integrate_shear
  
  curi = i
  curj = j
  curl = l
  myelg_calculate_shear = 9 *(CH0)**(-3) * OmegaM**2 *integrate(func, 0.01d0, 100.d0, 100 ,0)

end function myelg_calculate_shear

function myelg_integrate_shear(z)

  double precision, intent(in) ::z

  double precision :: myelg_integrate_shear, P, k, H

  k = curl/myelg_xi(z)
  !write(*,*) k
  P = myelg_powerspectrum(z,k)
  H = myelg_H(z)

  myelg_integrate_shear = myelg_g(curi,z)*myelg_g(curj,z)*(1+z)**2/H * P


end function myelg_integrate_shear

function myelg_OmegaM(a)
  double precision, intent(in) :: a
  double precision :: myelg_omegaM, H ,z
  z = 1.d0 / a - 1.d0
  H = myelg_H( z )
  myelg_omegaM = omegaM /(a**3)/(H)**2

end function myelg_OmegaM

function myelg_D(a)
  implicit none
  double precision, intent(in) :: a
  double precision :: myelg_D, igral

  PROCEDURE(function_template), POINTER :: func

  func => myelg_integrate_D
  
  igral = integrate(func,0.0001d0, a, 10000 ,1)

  myelg_D = a*exp(igral)

end function myelg_D



function myelg_integrate_D(a) 
  double precision, intent(in) :: a
  double precision :: myelg_integrate_D

  myelg_integrate_D = (1.d0/a)*(myelg_OmegaM(a)**gamma - 1.d0)

end function myelg_integrate_D



subroutine myelg_calculate_deltaH()
  double precision :: n, omegaL
  n = ns -1
  omegaL = 1 - OmegaM


  deltaH = 1.9E-5*exp(-1.01*n)*OmegaM**(-0.8 -0.05*log(OmegaM))*(1 + 0.18 * n*OmegaL)

end subroutine myelg_calculate_deltaH



function myelg_transfer(x)

  double precision, intent(in) :: x
  double precision :: myelg_transfer
  
  myelg_transfer = (log (1 + 0.171*x)/ (0.171*x)) *  (1 + 0.284*x + (1.18*x)**2 + (0.399*x)**3 + (0.49*x)**4)**(-0.25)


end function myelg_transfer

function myelg_powerspectrum(z,k)

  double precision, intent(in) :: z,k

  double precision :: myelg_powerspectrum, T, pi, D, a

  pi = 3.14159265
 
  T = myelg_transfer(k/keq)

  a = 1.0/(1.0+z)
  D = myelg_D(a)
  
   myelg_powerspectrum = 2*pi**2 *deltaH**2*(ch0)**(ns + 3) *   (k/ ( h))**ns * T**2 * (D/D0)**2

end function myelg_powerspectrum







function myelg_g(i, z)
  integer, intent(in) :: i
  double precision, intent(in) :: z
  double precision :: myelg_g, c
  PROCEDURE(function_template), POINTER :: func


!  N = 100

  xi = myelg_xi(z)

  func => myelg_integrate_g

  myelg_g = integrate(func, z_bins(i), z_bins(i+1), 1000 ,0)

end function myelg_g






function myelg_integrate_xi(z) 
  double precision, intent(in) :: z
  double precision :: myelg_integrate_xi

  myelg_integrate_xi =  1.0/myelg_H(z)

end function myelg_integrate_xi


function myelg_integrate_g(z) 
  double precision, intent(in) :: z
  double precision :: myelg_integrate_g

  myelg_integrate_g = (1 - xi/myelg_xi(z))*myelg_phi(z)

end function myelg_integrate_g



function myelg_xi( z ) 
  implicit none

  double precision, intent(in) :: z
  double precision :: myelg_xi
  PROCEDURE(function_template), POINTER :: func


  func => myelg_integrate_xi

  myelg_xi = Ch0*integrate(func, 0.d0, z, 150, 0)
  
end function myelg_xi

function myelg_H(z) 
implicit none
  double precision, intent(in) :: z
  double precision :: myelg_H, h
  
  h = OmegaM*(1+z)**3 + (1-OmegaM) *exp ( -3 * wa*z/(1+z))*(1+z)**(3* ( 1+w0 + wa))
  myelg_H = sqrt (h)

end function myelg_H


function myelg_phi(z)
  double precision, intent(in) :: z
  
  double precision :: myelg_phi

  myelg_phi = 1.5/z0 * (z/z0)**2 * exp(- (z/z0)**1.5)

end function myelg_phi





subroutine myelg_test_powerspectrum()

  integer :: i
  double precision :: l, z , P, k
  z = 0.5

  open(unit=10,file='ps.dat',action='write',status='replace')
  do i=1, 1000
     k = float(i)*0.01
     P = myelg_powerspectrum(z,k)
     !write(10,*) log(k), log(P)
     write(10,*) k, P
  end do
  close(10)


end subroutine myelg_test_powerspectrum


function integrate(func, xfrom, xto, N, flag)
implicit none
  integer, intent(in) :: N, flag
  double precision, intent(in) :: xfrom,xto
  PROCEDURE(function_template), POINTER :: func
  integer :: i
  double precision :: a, b, dx, c, integrate

  dx = (xto-xfrom)/float(N+1)
  c = 0
  if (flag==1) open(unit=12,file='int.dat',action='write',status='replace')

  do i=0, N-1
     a = xfrom + dx*float(i)
     b = xfrom + dx*float(i+1)
     c = c + (b-a) * 0.5*( func(a) + func(b)) 
     if (flag==1)write(12,*) a, func(a), c
  end do
  if (flag==1) close(12)
  integrate = c

end function integrate



      subroutine elg_testgiz(i, z)
        implicit none
        double precision :: val, om, w0, wa
        integer, intent(in) :: i
        double precision, intent(in) :: z
        integer          nbins
        parameter        (nbins=6)
        double precision zbins(nbins)
        external     ::    giz
        
        
        zbins(1)=0.01d0         
        zbins(2)=0.50d0
        zbins(3)=0.71d0
        zbins(4)=0.91d0
        zbins(5)=1.19d0 
        zbins(6)=2.5d0
        
        om = 0.3d0
        w0 = -1.0d0
        wa = 0.01d0
        call giz(om, w0, wa,zbins,nbins,i,z)
        !val=giz(om,w0dum,wadum,zbins,nbins,ibin,z(i)) ! linseeffektivitet, bin i 
        write(*,*) "G: ", val
        
      end subroutine elg_testgiz



subroutine test_lingrow

  double precision :: z, D1z, D10, a

  z = 2.00d0
  call LinGrow(z,OmegaM,w0,wa,D1z,D10)

  write(*,*) "ELG (D, D0) = ", D1Z, D10

  a = 1.d0/(1.d0 +z)

  write(*,*) "MINE (D, D0) = ", myelg_D(a), myelg_D(1.d0)


end subroutine test_lingrow



end module mod_myelg
